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ACCURACY  OF  THREE  UNCONDITIONAL!, Y-STABLE  FDTD  SCHEMES 
FOR  SOLVING  MAXWELL’S  EQUATIONS 


Guilin  Sun  and  Christopher  W.  Trueman 

Electromagnetic  Compatibility  Laboratory,  Concordia  University, 
7141  Sherbrooke  Street  West,  Montreal,  Quebec,  Canada  H4B  1R6 


Abstract  - This  paper  discusses  accuracy  limitations  due 
to  numerical  dispersion  and  time  step  size  for  three 
implicit  unconditionally-stable  FDTD  methods: 
Alternate-Direction-Implicit  (ADI),  Crank-Nicolson 
(CN)  and  Crank-Nicolson-Douglas-Gunn  (CNDG).  It  is 
shown  that  for  a uniform  mesh,  the  three  methods  have 
the  same  numerical  phase  velocity  along  the  axes,  but 
have  large  differences  along  the  diagonals.  The  ADI 
method  has  two  orders-of-magnitude  larger  anisotropy 
than  that  of  CN  and  CNDG.  CNDG  has  no  anisotropy 
at  certain  Courant  numbers  and  mesh  densities.  At  the 
limit  of  zero  spatial  mesh  size,  the  three  methods  have 
different  “intrinsic  temporal  dispersion”  for  a given 
time  step  size:  CN  has  no  anisotropy;  ADI  has  positive 
anisotropy  and  CNDG  has  negative  anisotropy,  which  is 
much  smaller  than  ADI.  The  Nyquist  sampling 
theorem  provides  a fundamental  upper  bound  on  the 
time  step  size  for  all  three  methods.  It  is  shown  that  for 
ADI  and  CN  the  practical  upper  bound  is  close  to  the 
Nyquist  limit,  but  for  CNDG  is  half  the  Nyquist  limit. 

Keywords  - Maxwell’s  Equations,  Finite-difference  time- 
domain  method,  Alternate-Direction-Implicit  method, 
Crank-Nicolson  method,  Douglas-Gunn  method, 
numerical  dispersion,  numerical  anisotropy,  accuracy, 
Nyquist  criterion. 


I.  Introduction 

Yee’s  Finite-Difference  Time-Domain  (FDTD) 
method  is  popular  for  solving  Maxwell’s  Equations  [1].  It 
is  second-order  accurate  in  both  time  and  space  [2].  Since  it 
is  an  explicit  method  [3],  it  is  easy  to  program  and  efficient 
to  run.  However  it  suffers  from  the  Courant-Friedrich-Levy 
(CFL)  limit  or  the  Courant  limit  on  the  time  step  size 
required  for  stability.  For  objects  with  fine  geometrical 
features,  using  a fine  mesh  size  greatly  reduces  the 
allowable  time  step  size,  which  causes  the  CPU  time  to  be 
prohibitively  long.  To  eliminate  the  CFL  limit, 
unconditionally-stable  methods  working  with  large  Courant 
numbers  are  desirable.  Early  in  1984,  Holland  [4]  proposed 
an  implicit  method  but  it  was  not  completely  stable.  In 
1995,  Shang  [5]  developed  an  efficient  characteristic-based 
algorithm  and  Fijany  [6]  proposed  a parallel  Crank- 
Nicolson  (CN)  method  by  decomposition  of  the 
eigenvalue/eigenvector  for  the  wave  equations  of  the 
second  order.  In  1999  and  2000,  Namiki  [7]  and  Zhen  et  al. 


[8]  suggested  an  Altemate-Direction-Implicit  (ADI) 
method.  In  2001,  Beggs  and  Briley  [9]  reported  a two- 
factor  scheme  by  combining  a characteristic-based 
approach  to  spatial  differencing  with  an  implicit  lower- 
upper  approximate  factorization  that  avoids  the  solution  of 
a tri diagonal  system.  Very  recently,  Sun  and  Trueman  [10] 
proposed  a Crank-Nicolson  scheme  with  Douglas-Gunn 
algorithm  (CNDG). 

This  paper  compares  the  numerical  dispersion  and 
anisotropy  of  the  ADI,  CN  and  CNDG  methods  used  in 
FDTD,  and  analyzes  some  of  their  common  characteristics 
and  differences.  For  simplicity,  Yee’s  mesh  [1]  in  2D  is 
used  with  a TEZ  wave  in  a linear,  isotropic,  non-dispersive 
and  lossless  medium.  This  paper  is  organized  as  follows. 
In  Section  II,  the  amplification  factors  are  listed  for  the 
three  methods.  In  Section  III,  their  numerical  dispersion 
relations  are  given  and  compared  for  a given  mesh  density. 
In  Section  IV  their  numerical  anisotropy  is  analyzed  and 
compared  at  different  Courant  numbers.  In  Section  V the 
time  step  size  limits  are  discussed  and  related  to  the 
Nyquist  criterion.  In  Section  VI  the  accuracy  limit  due  to 
dispersion  is  analyzed  for  zero  spatial  mesh  size. 

II.  Amplification  Factors 

The  update  equations  for  the  ADI-FDTD,  CN- 
FDTD  and  CNDG-FDTD  methods  are  listed  in 
Appendices.  The  ADI-FDTD  uses  two  sub-steps  [7] 
(Appendix  I).  The  first  sub-step  advances  time  from  step  n 
to  step  n+ 1/2.  The  field  component  £"+1/2  is  fully  implicit 

and  requires  solving  a tridiagonal  matrix.  The  second  sub- 
step advances  time  step  from  n+l/2  to  w+1,  and  the  field 

component  E”+]  is  fully  implicit  and  also  requires  solving 
a similar  tridiagonal  matrix.  The  time  step  n+ 1/2  is 
intermediate  and  the  field  values  at  this  step  are  non- 
physical. The  amplification  factors  for  the  two  sub-steps 
are  [7,11] 


l+r 


l + ri 


^-^tan-'ftl+rixi+r,2)-!) 


(1) 


£2  = 


1 + rx  -j  tan'1  «l+r,2)(l+r’)-l) 
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where  y = V-T  , rx=cAt  sin(PxAx  / 2)  / Ax  , ry- 

cAtsm{pyAy  12)1  Ay  , A/ is  the  time  step  size,  c~\f*jji£ 
is  the  physical  velocity,  s and  /j  are  the  permittivity  and 
permeability  of  the  material  respectively,  Ax  and  Ay  are 
the  spatial  meshing  sizes  along  x and  y axes; 
PX  = P cos(^),  py  = psm(</>),  p2=px+py  where  p is 
the  numerical  wave  phase  constant;  and  angle  (f>  is  the 
direction  of  travel  with  respect  to  the  x axis. 

The  CN  scheme  averages  the  field  components  at 
time  step  n and  n+l  to  maintain  second  order  accuracy  in 
time  [6,  10]  (Appendix  II).  It  has  no  intermediate  time  step. 
But  the  resulting  block  tridiagonal  matrix  is  very  expensive 
to  solve  by  direct  methods  such  as  the  Gaussian 
elimination  or  the  banded  matrix  method,  as  well  as  by 
iterative  methods  such  as  Successive  Over-Relaxation 
(SOR)  or  the  iterative  Alternate-Direction-Implicit  (IADI) 
method  [3].  Since  the  discretization  of  the  Poisson 
Equation  and  heat  equation  also  leads  tosu  ch  a block 
tridiagonal  matrix,  several  other  methods  to  solve  it  can  be 
found  in  Refs.  [3]  and  [12].  Finjany  et  al.  [6]  solve  the 
block  tridiagonal  matrix  by  eigenvalue/eigenvector 
decomposition. 

Different  solution  methods  for  the  block 
tridiagonal  matrix  may  have  different  amplification  factors 
and  different  numerical  dispersion  relations.  This  paper 
assumes  a direct  solution  method  of  the  block  tridiagonal 
matrix  such  as  the  Gaussian  elimination  for  the  CN-FDTD 
method.  The  amplification  factor  for  this  CN  method  is 
[10] 


± j tan' 


1 -ri-r* 


oCN 


(3) 


Note  that  the  magnitude  of  the  amplification  factor  for  the 
CN  method  is  unity. 

The  CNDG  method  avoids  an  expensive  direct 
solution  (Appendix  III).  It  factorizes  the  block  tridiagonal 
matrix  and  has  two  sub-steps:  the  first  sub-step  finds  the 
intermediate  field  value  H* , and  the  second  sub-step  gets 
the  solution  at  the  time  step  n+l  [10].  Thus,  like  the  ADI 
method,  the  CNDG  method  needs  to  solve  a tridiagonal 
matrix  at  each  sub-step.  The  amplification  factor  is 


± j tan' 


(1-r/XI-r/) 


oCNDG 


(4) 


magnitude.  On  the  other  hand,  in  ADI,  the  intermediate 
time  is  (n  + \/  2)At , and  each  sub-step  has  its  own 
amplification  factor.  The  magnitudes  of  the  two 
amplification  factors  for  the  two  individual  sub-steps  are 
reciprocals.  Their  product  is  the  magnitude  of  the 
amplification  factor  for  one  full  update  cycle,  and  is  unity. 
Note  that  the  unity  magnitude  of  the  amplification  factors 
ensures  that  the  three  implicit  methods  are  unconditionally 
stable  and  strictly  non-dissipative. 


III.  Numerical  Dispersion 


The  numerical  dispersion  relation  can  be  written 
for  the  three  methods  [10,  1 1]  as 

tan2  (coAt  / 2)  = rx  + r2  + grxr2 

7 7 / r\ 


where  the  factor  g is 


g = 


0 

- tan2  (coAt  / 2) 


ADI 

CN 

CNDG 


(6) 


Given  a frequency  co  , a mesh  size  Ax  and  Ay  , a time  step 
size  At , and  direction  of  travel  (f> , Eqn.  (5)  is  solved 
implicitly  for  the  phase  constant  p . For  a wave  to 
propagate  p must  be  real.  The  numerical  dispersion  is 
quantified  by  the  relative  velocity  that  is  defined  as  the 
ratio  of  the  numerical  velocity  to  the  physical  velocity. 


Fig.  1 Numerical  dispersion  with  mesh  density  100  and 
s=l,  5 and  1 0,  for  ADI,  CN  and  CNDG. 


Though  ADI  and  CNDG  both  have  two  sub-steps,  there  are 
some  differences.  In  CNDG,  the  time  level  of  the 
intermediate  step  is  unknown,  and  there  is  only  one 
amplification  factor  for  the  two  sub-steps,  of  unity 


The  numerical  dispersion  relation  in  Eqn.  (5)  for 
the  ADI  and  CNDG  methods  has  been  validated  with 
numerical  experiments  [10,11].  Fig.l  shows  the  relative 
velocity  as  a function  of  the  direction  of  travel,  using 
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Ax  = Ay  with  mesh  density  N = 100  cells  per  wavelength, 

at  different  Courant  numbers  s = cAt(  Ax  equal  to  1,  5 and 
10.  Note  that  the  numerical  dispersion  curves  of  the  CN 
method  and  the  CNDG  method  are  almost  identical,  and 
are  difficult  to  distinguish  visually  in  Fig.l.  The  relative 
velocity  along  the  axes  is  exactly  the  same  for  the  three 
methods,  but  along  the  diagonal  ADI  is  quite  different 
from  CN  or  CNDG.  Thus  their  numerical  anisotropies  are 
significantly  different,  which  will  be  discussed  in  the  next 
section. 

The  “grid-related  numerical  dispersion”  is  defined 
as  the  dispersion  at  zero  timest  ep  size.  Since  the  cross 
term  in  Eqn.  (6)  goes  to  zero  as  time  step  size 

approaches  zero,  the  grid-related  numerical  dispersion  is 
the  same  for  the  three  methods,  and  is  the  same  as  that  of 
Yee’s  FDTD  [2,  13]. 

IV.  Numerical  Anisotropy 

The  velocity-anisotropy  error  is  often  used  to 
quantify  the  numerical  anisotropy  [2],  In  an  isotropic 
medium,  the  wavefront  of  a cylindrical  wave  is  a circle, 
that  is,  the  phase  velocity  is  the  same  in  all  directions. 
However,  in  the  numerical  domain,  the  numerical  wave 
velocity  usually  depends  on  the  direction  of  travel.  Like 
Yee’s  FDTD  with  a uniform  mesh,  the  ADI  and  CN 
methods  have  the  largest  velocity  along  the  diagonals  and 
slowest  along  the  axes,  as  shown  in  Fig.  1. 

However,  the  CNDG  behaves  differently  from 
CN  or  ADI.  At  certain  combinations  of  the  mesh  and  time 
step  sizes,  the  velocity  u45 «,  along  the  diagonals  can  be 
slower  than  the  velocity  uQ « along  the  axes.  To  evaluate 
this  phenomenon,  the  following  definition  of  the  velocity 
anisotropy  is  used  for  a uniform  mesh 

A u = ■ x 1 00%  m 

min{w45o,u0o}  ^ 

This  error  definition  has  the  same  magnitude  as  that 
Taflove  and  Hagness  [2],  but  can  be  positive  or  negative, 
depending  on  which  velocity  is  larger.  Fig.  2 and  Fig.  3 
show  the  numerical  anisotropy  at  mesh  densities  from  50  to 
100  cells  per  wavelength  for  the  ADI,  CN  and  CNDG 
methods.  The  anisotropy  of  ADI  is  about  two  orders  of 
magnitude  larger  than  that  of  CN  and  CNDG.  As  the 
Courant  number  increases,  the  anisotropy  in  ADI  increases 
significantly.  The  CN  and  CNDG  have  the  same  behaviour. 
For  the  Courant  number  smaller  than  about  11,  the 
anisotropy  of  CNDG  is  always  smaller  than  that  of  CN;  for 
the  Courant  number  larger  than  11,  the  anisotropy  of 
CNDG  is  larger  than  that  of  CN  at  coarse  mesh  density  (for 
example,  50),  but  it  quickly  becomes  smaller  than  that  of 
CN  after  certain  mesh  densities.  For  instance,  if  the  Courant 
number  is  12,  at  the  mesh  density  50,  the  numerical 
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anisotropy  is  -0.148*0.01  for  CNDG  and  0.082*0.01  for 
CN;  at  mesh  density  55,  their  absolute  anisotropy  values 
are  the  same,  about  0.066*0.01.  However,  at  mesh  density 
60,  CN  has  0.055*0.01  but  CNDG  has  only  0.038*0.01. 


Fig.  2 Numerical  anisotropy  of  the  ADI  method  at  Courant 
numbers  1,5  and  10. 


Fig.  3 Numerical  anisotropy  for  CN  and  CNDG. 

Notice  that  the  anisotropy  changes  its  sign  for 
CNDG  at  certain  combinations  of  the  mesh  density  and 
Courant  number.  This  suggests  that  at  certain  Courant 
numbers,  a mesh  may  have  no  anisotropy.  Fig.  4 shows  the 
relation  between  the  Courant  number  and  mesh  density 
with  zero  anisotropy.  Note  that  neither  ADI  nor  CN  has 
this  behavior. 

V.  Time  Step  Size  Limitations 

In  the  numerical  dispersion  relation  of  Eqn.  (5),  if 
the  tangent  is  infinite,  the  phase  constant  becomes  complex 
for  all  the  three  methods.  This  limit  is  reached  when 


(8) 
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Fig,  4 The  relation  of  the  Courant  number  and  mesh 
density  at  zero  anisotropy  for  CNDG. 


Fig.  5 Time  step  size  limits  with  respect  to  mesh  density  for 
the  three  methods. 

Eqn.  (8)  is  recognized  as  the  Nyquist  limit  for  the  time 
sampling.  For  a given  mesh  density  N , the  Nyquist  limit 
can  be  written  in  relation  to  the  Courant  number  as 


N 


However,  Eqn.  (5)  implies  a smaller  time  step  size  limit 
than  the  Nyquist  limit.  For  example,  along  the  .x-axis, 
sin(/?Ax/2)  must  be  smaller  than  or  equal  to  one  if  p is  to 
be  real. For  the  ADI  and  CN  methods,  the  minimum 
velocity  always  occurs  along  the  axes,  so  the  limiting 
Courant  number  for  a give  mesh  density  can  be  found  from 
Eqn.  (5)  to  satisfy  tan(^y/ N)-s  . For  the  CNDG  method, 
at  larger  Courant  numbers,  the  minimum  velocity  is  along 
the  diagonals.  For  Ax  = Ay,  Eqn.  (5)  for  the  CNDG 
method  leads  to 


Fig.  6 Intrinsic  temporal  anisotropy  with  zero  mesh  size 
for  ADI. 


Fig.  7 Intrinsic  temporal  anisotropy  with  zero  mesh  size  for 
CNDG. 


, i-  1 + Jl~tan4(tt)A?/2) 

sin 2 (V2/7  . Ar/4)  = — * 

(s  tan(a)At  / 2)y  (10) 

where  P ^ is  the  phase  constant  along  the  diagonals.  If 
tan(<yA?  / 2)  > 1 , then  the  phase  constant  p becomes  a 
complex  number.  Thus  tan(&>AJ  12)  - 1 is  the  limit  for  the 
CNDG  method.  This  corresponds  to  s = N / 4 . The  three 
curves  shown  in  Fig.  5 are  the  Nyquist  limit,  and  the  step 
size  limits  for  the  ADI  and  CN  methods  and  for  the  CNDG 
method.  If  a method  uses  a Courant  number  larger  than  its 
limit  shown  in  Fig.  5,  numerical  attenuation  will  occur, 
which  does  not  correspond  to  the  physical  reality. 

VI.  Intrinsic  Temporal  Dispersion 

In  Yee’s  FDTD,  as  the  mesh  density  increases,  the 
numerical  dispersion  decreases.  In  the  limit  of  an 
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infinitely-fine  mesh,  that  is,  zero  spatial  mesh  size,  there  is 
no  numerical  dispersion  because  the  Courant  limit  forces 
the  time  step  size  to  be  zero.  Hence,  Yee’s  method 
collapses  to  the  continuous  case,  that  is,  no  discretization 
for  time  and  space.  But  for  the  unconditionally-stable 
implicit  methods,  because  there  is  no  Courant  limit,  fine 
mesh  size  can  be  accompanied  by  very  large  time  step  size 
and  the  method  remains  stable.  Even  in  the  limit  of  zero 
spatial  mesh  size,  the  methods  are  still  stable  for  any  time 
step  size,  but  there  is  numerical  dispersion.  The  numerical 
dispersion  at  zero  mesh  size  is  an  intrinsic  limitation  and 
may  be  termed  “intrinsic  temporal  dispersion”,  previously 
called  “time-splitting-related  dispersion”  [11].  The 
intrinsic  temporal  dispersion  is  different  for  the  three 
methods.  With  some  manipulation,  it  can  be  written  as 

u _ coAtll 
c tan(<yA/  / 2) 


method  is  two  orders-of-magnitude  larger  than  that  of  CN 
and  CNDG.  In  the  limit  of  zero  mesh  size,  ADI  and  CNDG 
have  anisotropy  but  the  CN  method  does  not.  Different 
from  the  ADI  and  CN  methods,  the  CNDG  method  may 
have  slower  velocity  along  the  diagonals  than  along  the 
axes,  and  may  have  zero  anisotropy  at  certain  combinations 
of  Courant  number  and  mesh  density.  The  three  methods 
have  time  step  size  limits  that  are  smaller  than  the  Nyquist 
criterion.  CNDG  has  the  smallest  time  step  limit.  The 
intrinsic  temporal  dispersion  is  a fundamental  accuracy 
limitfor  the  three  methods  and  is  much  larger  for  ADI  than 
for  CN  or  CNDG  along  the  diagonals. 

APPENDIX  I Update  Equations  for  ADI-FDTD 

The  ADI-FDTD  has  two  sub-steps.  The  first  sub- 
step is  advancing  time  from  step  n to  step  «+ 1/2  by  use  of 
the  following  update  equations  [7] 


Numerical  calculations  using  Eqn.  (11)  show  that  the 
intrinsic  temporal  dispersion  as  a function  of  the  direction 
of  travel  is  similar  to  that  shown  in  Fig.  1 but  smaller  in 
magnitude.  From  Eqn.  (1 1)  it  can  be  seen  that  the  relative 
velocity  is  not  a function  of  direction  of  travel  for  the  CN 
method.  Therefore  CN’s  anisotropy  is  zero  at  the  zero 
mesh  size  limit.  But  for  ADI  and  CNDG,  there  is 
anisotropy,  as  shown  in  Fig.  6 and  Fig.  7.  This  anisotropy 
is  termed  the  “intrinsic  temporal  anisotropy.”  Note  that  the 
ADI’s  anisotropy  is  about  30  times  larger  than  that  of 
CNDG  at  the  time  step  size  of  one-tenth  Nyquist  time  step 
size  limit. 

Note  that  in  Eqn.  (1 1)  the  tangent  cannot  be  larger 
than  unity  for  CNDG  method;  otherwise  the  velocity  will 
be  a complex  number.  This  time  step  size  limit 
corresponds  to 


This  is  half  of  the  Nyquist  limit,  and  coincides  with  that  in 
Section  V for  the  time  step  size  limit  of  the  CNDG 
method.  For  non-zero  mesh  size,  the  numerical  velocity  is 
always  smaller  than  the  intrinsic  temporal  dispersion, 
therefore  Eqn.  (1 1)  is  a fundamental  accuracy  limit  for  the 
three  methods. 


VII.  Conclusion 


This  paper  has  discussed  several  aspects  of  the 
ADI,  CN  and  CNDG  methods.  The  magnitude  of  the 
overall  amplification  factor  for  all  three  methods  is  unity; 
hence  they  are  all  unconditionally  stable.  The  numerical 
dispersion  is  the  same  along  the  axes  for  the  three  methods, 
but  differs  along  the  diagonals.  The  anisotropy  in  the  ADI 


Enx+V\i  + 1/2  J)  = E"(i  + 1/2  J)  + a,{ 

Hnz(i  + l/2,j  + \/2)-H"z(i  + l/2J-l/2)}/Ay  (1-1) 

Eny+U2(i,j  + 1/2)  = Eny(i,j +1/2)  -a,{ 

//”+1/2(/ +1/2, _/'  + 1/2) -//"+1/2(/- l/2,y'  + 1/2)}/ Ac  0-2) 

H"+U2(i  + 1/2,7  + 1/2)  = //”(/  + l/2,y  + 1/2)  + a2{ 

Enx(i  + \/2J  + \)-Enx(i+V2,j)}/Ay 
-a2{Eny*u2(i  + \J + 1/2) -En/U2(i,  7 + 1/2)}/ Ac 

(1-3) 

where  ax  - At  / 2s , a2  = At  / 2fu , At  is  the  time  step  size; 
e and/i  are  the  permittivity  and  permeability  of  the 
material  respectively;  Ax  and  Ay  are  the  spatial  meshing 
sizes  along  x and  y axes;  i and  j are  the  integer-number 
indices  of  the  computational  cells;  and  n is  the  time  step 
index.  In  this  step,  E”+u  2 is  implicit  and  can  be  found  by 
solving  a tridiagonal  matrix  of  the  form 

En/m  (i  - 1,7  + 1 / 2)  - {(Ac2 /ata2  + 2 )E"y*m  ( i , 7 + 1/2) 

+ E"/'n (f  + 1, 7 + 1 / 2)  = -(Ac2  / a, a2 )Eny  (f,  j + 1 / 2)  + 

(Ac/  a2  ){H"  (i  + 1 /2, 7 + 1 / 2)  — H"  (i  - 1 / 2, 7 + 1 / 2)} 

+ Ac/ Ay{E”  (i  + l/2, 7 + 1)  -E”(i  + 1/2,7)  + 

E”  (('  - 1 / 2, 7)  -E"(i-ll  2, 7 + 1)} 

(1-4) 

The  second  sub-step  advances  time  step  from  n+1/2  to  n+\ 
and  the  update  equations  are 
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1 ( / + 1 / 2,  j)  = 1 /2(/  + 1/2, J)  + fl,  { 

//;+1  (/  + 1 / 2,  j + 1 / 2)  - /-/:+l  (i  + 1 / 2,  j - 1 / 2) } / Ay  (j_5) 

Ey+i(‘J  + 1/2)  = Ey+lll(i,j  + l/2)-a,{ 

//;+l/2(i  + l/2,y  + l/2)-//;+1/2(i-l/2,y  + l/2)}/Ax  (1-6) 

H"+i(i  + 1/2,7  + 1/2)  = H”+ln(i  + 1/2,7  + 1/2)  + a2{ 

E?\i  + 1/2,7  + 1)  - ^r'O'  + 1/2,7)}/ Ay 
- flj{£j+1(l  + 1,7  + 1/2)  - £;+1(/,7  + 1/2)}/ Ax  (1-7) 

Since  £"+l  is  implicit  and  it  is  found  by  solving  a 
tridiagonal  matrix  of  the  form 

E^a  + 1/2,7  - l)-{Ay2/fl,a2  + 2)  Era  + 1/2,7)  + 

Enr  a + 1/2,)  + 1)  = -(Av2  / a,fl2  )£f' /2(*  +1/2,  J)  + 

(Ay/a2){//;+,/2(/  + l/2,7-l/2)- 
/f"+,/2(j  + 1/2,7  + 1/2)} 

+ Ay/Ax{£'"+1/2(/  + l,7+l/2)-£'"+,/2(i,7  + l/2)  + 

£;+,/2(< +i,y  - 1/2)  - Enr'\uj  - 1/2)}  (I.8) 

APPENDIX  II  Update  Equations  for  CN-FDTD 

The  CN  scheme  averages  the  field  components  at 
time  step  n and  n+ 1 to  maintain  second  order  accuracy  in 
time  as  follows  [6,  10] 

Erx+\i+V2,j)=J?x(.i+VZj)+ 

a,  {I?*'  (/ + 1 / 2, 7 + 1 / 2)  - /^+1  (i  + 1 / 2, 7'  - 1 / 2) + 

//;(/+l/2,7+l/2)-/^(/+l/2,7-l/2)}/Ay  (nI) 

£;+i(/,7'+1/2)  = £;(/,  7+1/2)- 

//;+1(/+1/2,7'+1/2)-//;+i(i-1/2,7+1/2)+ 

H"  (i'+l/2, 7+l/2)-//;(/-l/2,7+l/2)}/Ax  (n*2) 

//;+l  (/  + 1 / 2, 7 + 1 / 2)  = H"  (i  + 1 / 2, 7 + 1 / 2)  + 
a2  {£;+1  (i  + 1 / 2, 7 + 1)  - £;+1  (i  + 1 / 2, 7)  + 

O'  + 1/2,7  + 1)  - O'  + l/2,7)}/Ay  - 
a2{£;+1  a + 1, 7 + 1 / 2)  - £;+l  (i,7  + 1 / 2)  + 

5;  (/  + 1,  j + 1 / 2)  - e;  (i,  / + 1 / 2)}  / ax  (n_3) 

The  block  tridiagonal  matrix  to  be  solved  for  the 
CN  method  is 


/^+1(/+l/27'+l/2)-fl,a2{/^+l(i+3/2,7+l/2)+ 
/^h'i(i-1/2,7+1/2)-2/^,+1(i +1/2,7 +1/2}/ Ax2 
-al£j2{/^+l(/+l/27+3/2)+/^+l(/+l/2,7-l/2)- 
2HTa +1/2,7  +1/2)}/ Ay2  = 
O+l/2,7+l/2)+a102{/^(i+3/2,7+l/2)+ 
0-l/2,7+l/2)-2/^(i+l/2,7+l/2)}/Ax2 
+ala2{f^(i+\/Zj+3l2)+H”2(i+l/2,j-\/2)- 
2/^(i +1/2,7 +1/2)}/ Ay2 

+2a2{K(1’+l/2,7+l)-^x(< +1/2,7)/ Ay-]- 

t^(/+l,7+l/2)-£j(/, 7+1/2)]/ Ax} 


APPENDIX  III  Update  Equations  for  CNDG-FDTD 

The  update  equations  of  CNDG-FDTD  are  the 
same  as  the  CN  method.  The  difference  between  the  CN- 
FDTD  and  CNDG-FDTD  is  the  method  used  to  solve  the 
block  tridiagonal  matrix.  In  CNDG-FDTD,  it  is  factorized 
into  two  sub-steps  [10] 

//*  (/  + 1 / 2,  y + 1 / 2)  -ajfl2  {H*z  (z  +3  / 2,  y + 1 / 2) + 

//*(/-l/ 2,7 +l/2)-2//^(i +1/2, y'+l/2)}/ Ax2 

=H"(i+l/2,j+l/2)+ala2{H”(i+3/2,j+l/2)+ 

//"(/-l/2,7+l/2)-2//”(/+l/2,7+l/2)}/Ax2 

+o,a2  {H"  (/+1/2, 7+3/2)+//;  (i+l/2, 7-1/2)- 

2//;(/+l/2,7+l/2)}/Ay2 

+ 2a2  { [£"  (/ + 1 / 2, 7 + 1)  - £"  (/ + 1 / 2, 7)]  / Ay  - 

t£;(/+l,7+l/2)-£:;(/,7+l/2)]/Ax}  (in-1) 


//;+l(i  + 1/2,7  + 1/2)  - + 1/2,7  + 3/2)  + 

//f  (i  + 1 / 2, 7 - 1 / 2)  - 2//;+1  (/  + 1/2, 7 + 1/2)}/  Ay2 
= //*(/  + 1 / 2,7  + 1 / 2)  (III-2) 

At  each  sub-step  a tridiagonal  matrix  must  be 
solved.  It  can  be  seen  that  the  ADI  method  solves  the 
electric  fields  implicitly,  but  the  CN  and  CNDG  methods 
solve  the  magnetic  field  implicitly  in  the  2D  TEZ cas  e. 
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